/*Calculo de la energia y fuerza dipolar de forma directa*/
#include <stdio.h>
#include <math.h>
#include <fstream>
#include <iostream>
#include <string.h>
#include <iomanip>
#include <stdlib.h>

using namespace std;
#define NCUT              5    /* number of replicas */
#define NP                 8    /* number of particles: */
#define L                 10.0  /* system length */
#define COULOMB_PREFACTOR  1.0  /* bjerrum_length*temperature */

/*Definimos los arrays de posicion de particulas y fuerza*/
double x[NP], y[NP], z[NP], Q[NP], force[NP][3];

/* Function Headers */
double directsum(void);

int main() {
	double s1 = 0.50*L;
	double E_Coulomb;     /* Coulomb-energy of the box */

	/* define the particle coordinates and charges such such that this
	gives an NaCl-lattice upon periodic replication: */

	x[ 0] = 0;  y[ 0] = 0;  z[ 0] = 0;  Q[ 0] =  1.0;
	x[ 1] = s1; y[ 1] = s1; z[ 1] = 0;  Q[ 1] =  1.0;
	x[ 2] = s1; y[ 2] = 0;  z[ 2] = 0;  Q[ 2] =  -1.0;
	x[ 3] = 0;  y[ 3] = s1; z[ 3] = s1; Q[ 3] =  1.0;
	x[ 4] = 0; y[ 4] = s1;  z[ 4] = 0;  Q[ 4] = -1.0;
	x[ 5] = s1;  y[ 5] = 0; z[ 5] = s1; Q[ 5] = 1.0;
	x[ 6] = 0;  y[ 6] = 0;  z[ 6] = s1; Q[ 6] = -1.0;
	x[ 7] = s1; y[ 7] = s1; z[ 7] = s1; Q[ 7] = -1.0;

	for(int i=0;i<8;i++)
		cout<<x[i]<<"  "<<y[i]<<"  "<<z[i]<<"  "<<Q[i]<<endl;

	for(int i=0;i<NP;i++){
		for(int k=0;k<3;k++){
			force[i][k]=0.0;
		}
	}

	E_Coulomb = directsum();
	fprintf(stderr, "Coulomb energy (Direct Sum) is %20.15le\n", E_Coulomb);

	for(int i=0;i<NP;i++)
		cout<<force[i][0]<<"  "<<force[i][1]<<"  "<<force[i][2]<<endl;

	return 0;
}

double directsum(void) {
	int i, j;
	/**********************
	** check system
	*********************/
	double qsum = 0.0;
	for (i = 0; i < NP; i++) {
		qsum += Q[i];
		if (x[i] < 0.0 || x[i] > L)
			printf("WARNING: X-coordinate of particle %d is not in the central image!\nx[%d]=%f\n",  i, i, x[i]);
		if (y[i] < 0.0 || y[i] > L)
			printf("WARNING: Y-coordinate of particle %d is not in the central image!\ny[%d]=%f\n",  i, i, y[i]);
		if (z[i] < 0.0 || z[i] > L)
			printf("WARNING: Z-coordinate of particle %d is not in the central image!\nz[%d]=%f\n",  i, i, z[i]);
	}

	/**********************
	** compute sum
	*********************/
	double E;
	double qij;
	double nq;
	int nx,ny,nz;
	double r,r2,rx,ry,rz,rnx,rny,rnz;
	double Co;
	double NCUT2=NCUT*NCUT;
	double mad=0.;
	double coef;
	E=0.0;


	for(i=0; i<NP; i++){
		for(j=0; j<NP; j++){
			qij=Q[i]*Q[j];
			rx=x[i]-x[j];
			ry=y[i]-y[j];
			rz=z[i]-z[j];

			for(nx=-NCUT;nx<=NCUT;nx++){
				for(ny=-NCUT;ny<=NCUT;ny++){
					for(nz=-NCUT;nz<=NCUT;nz++){
						nq=nx*nx+ny*ny+nz*nz;
						if(nq < NCUT2){
							if((nq!=0) || (nq==0 &&i!=j))
							{
							rnx=rx+nx*L;
							rny=ry+ny*L;
							rnz=rz+nz*L;

							r2=rnx*rnx+rny*rny+rnz*rnz;
							r=sqrt(r2);
							Co=rnz/r;
							E += (3.-Co*Co)/(r2*r);
//							coef=qij/r2;
//							force[i][0]+=coef*rnx;
//							force[i][1]+=coef*rny;
//							force[i][2]+=coef*rnz;
							}
						}
					}
				}
			}
		}
	}

	E *= 0.5;

	return E;
}

